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SUMMARY 

This paper is concerned with modeling flexural and membrane type waves 
existing in various submerged (or in vacuo) plate and/or shell finite element 
models that are excited with steady state type, harmonic loadings proportioned 
to e^ u t. only thin walled plates and shells are treated wherein rotary inertia 
and shear correction factors are not included. More specifically, the issue of 
determining the shell or plate mesh size needed to represent the spatial dis- 
tribution of the plate or shell response is of prime importance towards suc- 
cessfully representing the solution to the problem at hand. To this end, a 
procedure is presented for establishing guide lines for determining the mesh 
size based on a simple test model that can be used for a variety of plate and 
shell configurations such as, (1) cylindrical shells with water loading, (2) 
cylindrical shells in vacuo, (3) plates with water loading, and (4) plates in 
vacuo. The procedure for doing these four cases is given, with specific numer- 
ical examples present only for the cylindrical shell case. 


INTRODUCTION 

This paper addresses the topic of modeling flexural waves and membrane 
waves present in various types of shell and plate type structural configura- 
tions. The issue at hand is arriving at a simple procedure for determining a 
mesh size adequate to represent the details of the spatial response distribu- 
tions necessary to achieve some desired level of accuracy. To be sure, it 
would be too large an undertaking to answer this question for all possible 
plate and shell configurations that may arise, however, the selected class of 
thin walled cylindrical shells and flat plates are often the major building 
blocks of a good deal of structures. Therefore, the paper will focus on these 
two configurations, with the major emphasis on the cylindrical shell employing 
CCONEAX elements with ax i symmetrical loading. Physical problems with finite 
length dimensions exhibit solution responses that often have the form of 
standing wave patterns as a result of reflections from the shell boundaries. 
Further, these solutions do not have a single clearly defined constant ampli- 
tude traveling wave component as one would have, in say, an infinitely long 
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shell or plate. A constant amplitude traveling wave (prop$aating in the x 
direction at frequency to) response, R, of the form R = R^ 1 v-Yx-hjt) a <j e sir- 
able form to seek because the constant amplitude Rq and associated spatial wave 
length ^ Y = 2 tt / y can then be used as a measure for determining the ability of a 
particular element to represent the desired wave propagation phenomenon. Devi- 
ations of the finite element amplitude response from the exact response con- 
stant R and/or deviatons of the finite element response phase angle period 
from the exact spatial period can be used to judge the mesh refinement neces- 
sary to achieve a desired level of accuracy. The thrust of this paper is not 
to give a hard and fast rule for the number of elements per wavelength neces- 
sary to achieve some desired level of accuracy, but rather to provide a proce- 
dure for allowing one to establish his own level of accuracy. We take this 
approach because rules of thumb are often dangerous, particularly in the area 
of wave propagation in plate and shell structures. There are cases, for exam- 
ple, in a plate or cylindrical shell where a cutoff point exists such that the 
particular problem parameters (geometry, frequency, and physical constants) 
result in a situation where there is no traveling wave. If the problem param- 
eters happen to be such that the traveling wave root is close to the cut-off 
point, a finer mesh size might be needed to properly represent the propagating 
wave situation. 

The symmetrically loaded (e independent loading) infinitely long cylinder 
shell (with or without fluid external fluid present) is selected as the model 
for examining flexural and membrane traveling waves (see Figure 1). This same 
model can be used to treat all four plate and shell cases (with and without 
water) discussed above. The plate cases can be realized by letting a -*■ 00 and 
the in vacuo cases achieved by setting the density of the fluid equal to zero. 


SYMBOLS 


a 

C 





D 

f 



radius of cylindrical shell (in) 

fluid sound speed (in/sec) 

shear velocity = G/p s ' psi 

in vacuo plate wave speed, n/e/(o s (1-v 2 )) (in/ sec) 

propagating wave phase velocity, y/y 

in vacuo plate wave speed, ^lhp^//D (in/ sec) 

plate modulus = Eh3/[12 (1 - v 2 )] of rigidity (lb/in) 

harmonic frequency (Hz) 

shear modulus of elasticity (psi) 

Hankel Bessel functions of the first kind ord;r 0 and 1 
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h = plate or shell thickness (in) 
i = nTT complex number 

K 0 ( ) ; K^O = modified Bessel functions of order 0 and 1 
k = w/c acoustic wave number (in - !) 
kfp = w/Cfp in vacuo plate wave number (in-1) 
kp = to/ Cp in vacuo plate wave number (in-1) 

L = finite element model length (in) 

M z = shell line moment (Ib/in/in) 

N z = axial line membrane force (lb/in) 
p = fluid pressure (psi) 

Q = transverse shell line shear force (Ib/in) 
r = radial cylindrical coordinate (in) 
t = time variable (sec) 
u = axial motion of plate or shell (in) 

U Q = amplitude of u (in) 
w = radial motion of plate or shell (in) 

W Q = amplitude of w (in) 

x = plate coordinate direction normal to plate (in) 
z = axial cylindrical coordinate variable (in) 

= wavelength of propagating wave 
Y = traveling wave number for z direction (in-1) 

Y r = real part of y 
Y- j = imaginary part of y 
p = mass density of fluid (lb/sec2 in-4) 

P $ = mass density of structure (lb/sec2 in _ 4) 
v = Poisson's ratio of structure 
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a) = 2nf = harmonic frequency (rad/sec) 
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ANALYTICAL WAVELENGTH FORMULATION FOR 
PLATES AND CYLINDRICAL SHELLS 


The analytical formulation employed for determining the exact relationship 
for freely propagating waves in plate and cylindrical shell structures is 
important for its own sake, however, the formulation also directly leads to the 
procedure employed to set up the NASTRAN wavelength accuracy demonstration 
models. Consequently, some of the important details of the development for the 
freely propagating wave characteristics equation are given. The procedure 
given will closely follow the one given in ref. 2, except that freely propagat- 
ing rathe" than standing wave configurations are considered by changing the 
axial variation as a complex exponential variation in z rather than a cosine 
variation. 


Of eventual interest is both the flexural wave and axial (membrane) wave- 
lengths for plates (with and without water present) and for infinitely long 
cylindrical shells (with and without water present). This constitutes 2 x 2 x 
2=8 different situations. However, we can treat the eight cases at the same 
time by first treating the most complex case of the submerged infinitely cylin- 
drical shell. By taking appropriate limits of this case (i.e., water density 
-*•0 or shell radius ■> ») we can recover the other cases of interest without 
requiring any new analyses. 


We start with the governing simplified thin wall shell field equations for 
the cylindrical shell (ref. 1) and corresponding wave equation for the external 
fluid, namely 



v_ 

3w 

0 2 u 

1 - n 

9^ + 

a 

3z 

~ w * 

C P 

v 9u 


w j. 

h 2 9"w 

. 9 2 W 

a 9z 

+ 

a 2+ 

17 9? r 

+ w 

V 2 p - 

1 9 2 

C T 9t 

£ = 0 



( 1-v 2 )/ (Eh) 
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( 2 ) 


subject to the usual boundary condition relation relating the pressure gradient 
normal to the shell and the shell acceleration. 


9p _ 9 2 w 

9r " p dV 


(3) 


r=a 

where u, w are the axial, and radial displacements of the shell (see Figure 2a 
for positive sense); p is the pressure in the fluid; h the shell thickness; a 
the radius of the cylindrical shell; r, z are the radial axial cylindrical 
coordinates; E,v are Young's modulus and Poiss on's ratio of the shell; Cp is 
the plate wave speed parameter ( n/T7Tp s (1 -v^) ) ); and c is the acoustic wave 
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speed of the fluid. The shell equations (1) are for classical thin wall plate 
and shell theory and do not have rotary inertia and shear correction factors 
included. Consequently, the area of interest focused on herein will be in a 
frequency range of w such that the above two methods correction factors are not 
necessary. This point will be expanded upon later in the paper. 


The 9 symmetrical shell motion representing propagating waves in the +z 
direction are assumed to take the form 

u = U e“ lwt 

o 

v = 0 (no hoop motion) (4) 

w = e iY2 e"^ 

where and W 0 are the yet to be determined wave amplitudes, and y is t'te 
propagating wave number. For outward going waves, the form of the pressure 
solution identically satisfying the wave equation (2) is given by ref. 2, 

p(r,z) = P Q H^(r *Jk 2 -y 2 1 )e^ YZ e -lti)t (5) 

which can be easily verified by direct substitution of equation (5) into equa- 
tion (2), where Po is the yet unknown pressure amplitude, and H^)( ) is the 
Hankel, function of the first kind of order zero. A similar expression for p 
with Hfl ( ) (Hankel functions of the second kind) also satisfies the wave equa- 
tion, but represents inward coming waves (when k2 >y2) 0 r results in an ever 
increasing exponential increasing pressure field with increasing r when k2 < y2„ 
Neither of these situations corresponds to the physical problem at hand, thus 
only the Hankel function of the first kind is retained. 


The characteristic equation resulting in an interaction equation relating 
the driving frequency wand admissible propagating wave numbers, Y, is obtained 
by substituting equations (4) and (5) into equation (1), subject to the bound- 
ary condition (3): actually substituting equations (4) and (5) into condition 

(3) leads to the relation between the surface motion amplitude W Q and the pres- 
sure amplitude P Q , namely 

P n = "P w n 0)2 (5a) 

N/k 2 -y 2 ' H ( j l) (a k 2 -y 2 ') 


Thus the resulting two linear equations in the W Q , U Q coefficients is given by 


UJ/Sr-Y‘1 +H 




V (f) + “o F * T2 Y " ‘ ( “ /C p )! 


( 6 ) 


(1-v 2 )poj 2 hW (a 

Eh *J k 2 -y 2 h} 1 ) (a k 2 -y 2 ') 


= 0 


The nontrivial solution is obtained by setting the determinant of the 
amplitude coefficients equal to zero resulting in 
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- r 


± - ( W /Cp ) 2 - 


for > y2 


Eh’sTM-y 2 ^ H^a *1 k 2 -y 2 ' ) 


(7a) 
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or in another equivalent form of 


CO 2 

rr - Y 2 


Jr + £ y- - Wc n |> - 


Ehs/r* 7 ^ Ma^-k 4 ' ) 




(7b) 


for y 2 > k 2 


where the Kq ( ) and Kj() are modified Bessel functions and are related to the 
Hankel functions by the identity 


K (X) - \ i n+1 H^(iX) 


Noting that the in vacuo flexural wave number, kp , formula for a plate is 


given by 


k f P - 


(8) 


(9) 


and the in vacuo membrane compressional wave number, kp , given by 

k = — 
p C P 


to 




do) 




and finally the acoustic wave number by 


k = “ 
c 


( 11 ) 


where c is the acoustic wave speed of the fluid medium; the characteristic 
equation for the traveling wave number, non-dimensional ized with respect to the 
acoustic wave number k, can algebraically be rewritten as 


((kp/k) 2 -(y/k) 2 ) 
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The original dimensional form of the characteristic equation (7) depended on 
the eight parameters h, w, E, v> pi p s> c, for the resulting wave numbery; 
however, the above equation (12) fomronly depends on the reduced number of 5 
non-dimensional parameters, namely lu/k , ka, p/p St a/h, and v for the result- 
ing non-dimensional wave numbers y/k. The term (kfp/k) is not an independent 
parameter since it can be expressed in terms of the other ones through the 
relation 


■4 E = (12)* ^(y 10(f)/ (kaV 


(13) 


Roots of Characteristic Equation and Mode Shapes 


The roots of the frequency equation (12) cover several special cases and 
require further elaboration. The roots of interest for freely propagating 
waves are those for which y = y~ is a real quantity, wherein the shell response 
has arguments (see equation (4 ) ) of the traveling wave form i ( y r z - wt) and 
corresponds to a wave propagating in the +z direction and the corresponding 
wavelength is aiven by 

* - ^ (14) 

Possibly roots of the characteristic equation are complex (e.g., this could be 
the case when employing equation (12a) due to the complex form of the Hankel 
functions H n ( ) = J n ( ) + i Y n ( )). In these situations, the complex root can 
be written in the form y = y r + iyj and substitution of this quantity into the 
sf. ?11 response .equation (4), shows that the axial z response is proportional to 
gHYr+iYiJz _ e iy r / g-y-jz . Thus the solution amplitude would either reduce 
(Yi > 0) or grow indefinitely (yj < 0) with z according to the sign of yi. The 
wavelength of the reducing (or growing) fluxuations is still given byx r = 2Tr/y r 0 


The case of interest in the remainder of this paper is the one for which 
there are freely propagating waves of constant amplitude. This situation will 
result when (Y/k)2 > 1 wherein the modified Bessel functions Kq and Kj have 
real arguments and consequently equation (12b) usually results in real roots 
for y. The situation (y/k)2 > l implies the acoustic wavelength X a = 27t/k is 
longer than the propagating wavelength, X r = 2 y r 

For cylindrical shells, unlike the infinite plate, the axial membrane 
force N z , will also result in a radial component of motion. Consequently, the 
axial and radial motions for the membrane propagation or the flexural wave 
propagation are coupled. In fact, the same characteristic equation (12b) can 
be used to determine both the membrane and flexural propagating wave numbers 
Yax> Yf l . Whether a given root.Y , corresponds to the membrane or flexural wave 
can be established by examining the mode shape. More specifically, the ratio 
of U 0 /W 0 can be solved from either of the two homogenous equations (6); upon 
substituting the root, y, into say the first of equations (6) yields 
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(15) 



- v(y/k) e lTT ^ 2 


Upon subsituting the y/k root in question into equation (15), the size of Uo 
relative to Wq can be established, wherein IUq/%! >> 1.0 implied a dominant 
axial motion hence signifying a membrane propagating wave and |Ub/%l << 1*0 
implies a dominant radial motion hence signifying a flexural propagating wave. 


Limiting Cases of General Characteristic Equation 

Several limiting cases of the general characteristic equation (12b) not 
only provide useful information by their own merit (i.e., the traveling wave 
characteristic equation of a plate in fluid or air; or, the case of a cylinder in 
air) but at the same time provide insight with regard to where (i.e., the range 
of (y/k)) to search for the root. It is instructive to start with the simplest 
case of an infinite flat plate in vacuo and build up to the more general case 
of a submerged infinite cylinder. 

• infinite plate in vacuo 


This situation is realized by taking the limit as the fluid to structure 
mass ratio goes to zero, p/p s -*■ 0, and as the non-dimensional frequency param- 
eters ka -*■ °° . It is noted that ka 00 can be realized by having the radius 

a ■> » at a finite frequency w = k/c. Thus equation (12b) reduces to 



(16) 


which by inspection has the roots 


y = k and y = k 


fp 


Substituting y = kp into equations (6) shows that 
Lb K0) + W 0 • (vkp/®) = 0 
U Q 1(0) + W Q • Ch2 kj/12_ j^(o»/ Cp )2] = 0 

~ To'""" 


(17) 


thus, W 0 = 0 and U 0 is any constant, and hence this mode corresponds to pure 
axial motion with no radial motion, thus indicating the y = k_ root is for the 
membrane traveling wave. Similarly, substituting y = kf into equations (6) 
yields 


t 0 


U Q i [u)2/ c 2 - k$ ] + W Q • (0) = 0 


V [0] + W 0 • [0] = 0 


(18) 
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thus t =0, and is any constant and hence corresponds to pure flexural 
motion with no axial motion. 

• infinite plate in fluid 

This situation is similarly obtained by passing to the limit as a -*■ 00 in 
equation (12b) and noting that (Kq ( x)/K. (x) ) -+ 1 as x -*-«>. Thus equation (12b) 
reduces to 



where again separate membrane traveling waves and flexural traveling waves result 
by separately equating the [ ]j and [ ] 2 terms to zero in quation (19). Thus 
by inspection, as in the in vacuo plate case, 

i’k- 


is again a root for the membrane propagating wave and the root 


1 



(p/ P S ) 

khN, (jc)* _1 



of 


(19a) 


provides the flexural wave root. It is noted that equation (19a) is algebrai- 
cally equivalent to equation (7.10) of ref. 2. Upon inspecting equation (19a), 
it is seen that the effect of the water presence is to make the Y/k root a 
larger value (i.e., the flexural wavelength is smaller) than what would have 
been the Ca:e in the absence of external water. Thus the water has no effect 
on the membrane wave root, however the flexural wave root is effected and must 
be solved numerically with some sort of root searching scheme. 


• infinite in vacuo cylindrical shell 


The characteristic equation is again obtained from the general (12b) case, 
by passing to the limit as (p/p s ) +0, thus obtaining 


lk P 


- Y 




2 2 

¥- = ° 


( 20 ) 


which corresponds to a cubic equation in the desired y2 root and consequently 
can be solved for exactly without resorting to i numerical root finder. Due to 
the presence of the last term in equation (20), y = kp is no longer the exact 
root for the membrane traveling wave, however, the (vy /a)2 does not usually 


shift the membrane wave number root very far from the plate solution, y - kpas 
can be seen after solving for the exact roots of equation (20). Only positive 
values for y 2 obtained from the cubic solution will correspond to freely propa- 
gating waves. Although the exact membrane and flexural roots are coupled 
through the last term in equation (20), an approximation of the root, say for 
roughing out the mesh size needed, can be obtained from the approximations 
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obtained after dropping the normally small (v y/a)2 term, hence 


for membrane wives, and 

y . k fp "TyF 

for flexural waves. 





MODELING FLEXURAL AND MEMBRANE WAVES WITH NASTRAN 

A simplified model Is needed for the determination of the necessary number 
of elements per wavelength required to accurately model the membrane and/or 
flexural freely propagating waves In NASTRAN type elements. The procedure used 
here is to cut a segment (say at least one wavelength long) out of the infi- 
nitely long model. At one cut of the model (call It the starting end "A"), 
moments and forces existing in the freely propagating wave are applied explic- 
itly; at the other end (call it the termination end "B"), an appropriate 
absorbing boundary conditions is applied (dampers that relate shell edge veloc- 
ities to edge moments and forces). The boundary absorbers Implicitly apply the 
appropriate amount of moment or force that would have been present internally 
at that section of the shell. In theory, it would be possible to apply explic- 
it appropriately phased mements and forces at the termination end as well, how- 
ever, this is not done for two reasons. First, there Is a certain amount of 
phase angle drift (an amount beyond the expected termination phase of yL) that 
exists as the wave propagates along the axis of the shell (or plate) and it Is 
not known a priori what the phase drift is so that the termination moments and 
forces can't be appropriately adjusted; the application of boundary absorbers, 
however, does not require any phase drift adjustment. Secondly, the implemen- 
tation of the absorbers gives future insight on how to truncate finite element 
tells or plates in such a manner that large problems requiring premature model 
termination can be made. 


In Vacuo Cylindrical Shell 


• determination of edge loading moments and forces 


The appropriate moments and forces at the starting end A can be citermlned 
from the relationships relating the shell motion u, w to the shell edge loads 
N^, QA, m^ and a'e given by the shell theory relations (ref. 1) 





v ? 


w) 

a / 


(21) 


where the equation (21) moments and forces are line loads (l.e., loads per unit 
circumferential arc length) obeying the shell sign convention of Figure 2a. 
Upon substituting equations (4) into equations (21), it follows that 
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mJ = - y 2 D W q e 1 (yz-wt) qA = _ i y 3 D Wq gi(YZ-u)t) 


( 22 ) 


where the relation between U 0 and Info is governed by equation (15) for cylindri- 
cal shells. For plates (a -*■ <») , Uq and W 0 are independent, where for membrane 
waves U 0 denotes the wave amplitude with ^ = 0 and for flexural waves Wo is 
the wave amplitude with J Q = 0. 


Equations (22) are not in a form ready to use in NASTRAN CCONEAX elements 
for three reasons: first, NASTRAN input is not in the usual line load form as 

is the case with shell theory, but rather are total values wherein the li-e 
loads must be multiplied by 2Tra; second, the finite elemenv type sign conven- 
tion (like quantities at both ends of the member are positive in the same 
sense) is different from the shell theory convention (e.g., compare Figures 2a 
and 2b) and; third, NASTRAN has an e +la)t hardwired into the code, consequently 
an adjustment in the analytical solution is necessary to compensate for this 
point. Specifically, a wave traveling in the +z direction also can be repre- 
sented with an i(-yz + ojt) type argument, thus by replacing y with -y and w 
with -a) in the analytical solutions will accomplish the corresponding correc- 
tion. Thus accounting for all three compensations and also incorporating 
equation (15), the NASTRAN loading for the cylindrical shell at the starting 
cut A, evaluated at z = 0, is given by: 


2! 
M 3=> 

II 

2-rra 

y z D 

W Q e +iajt 

NASTRAN 





cut-A 

q a = 

2-na 

y 3 D 

.. +90° i +icot 

W 0 e e 

loading 

for 





cylinder 


2ua 
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IW 




(23) 


where W Q is any suitable displacement amplitude factor that is selected by the 
user, and y f$ the root directly from the appropriate frequency equation 
(sign changes in y and w have already been compensated for in setting up equa- 
tion (23)). The e+iwt factors are of course omitted when entered data into 
NASTRAN since this factor is automatically accounted for internal to the 
program. 

• determination of shell edge termination absorbers 

The membrane or flexural wave incident upon the shell termination will 
reflect from the end unless the appropriate boundary condition is inserted to 
simulate the effect of having an infinitely long shell. Here the appropriate 
boundary absorber is developed for each of the three shell edge loads. Follow 
ing the same procedure used to setup equations (23), the internal moments and 
forces at the termination end are given by 
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nB _ o a 3 n n +i 90° u +iujt o -iYL 

Q = -2iTa y D e W Q e e 

~B -27ia h p c a) 2 e +,90 ° M _+iojt -iYL 
N = Y U Q e e 


(24) 


w * 


The equations (24) can be written in terms of velocity, by noting u = + iwu 
i w w, dw/dx = 6 = ( - i y) ( icow) thus equations (24) become 


S| . (6) 


rB -2TiaY 3 r ,• . 

q * — £- m 


(25) 


where the above forces are the forces on the shell edge B. Dampers attached to 
the end of the shell would have an internal force or moment equal and opposite 
to the value at the shell edge, thus the forces and moments internal to the 
damper ’re given by 


^ = C 
d e 


( 9 ) 


Q » C r 


where the damper values are given by 
2-rra D 


r - 2tra D y 2 
c y 

where the propagation phase velocity, 


V 


U) 


(w) W = • (u) (26) 


{ NA^TRAN 

CCONEAX (27; 

DAMPERS 

Cy, is defined by 

(28) 


where y is the root of the characteristic equation. In the case of a plate, 
the same dampers can be used, ex r 2 pt the 2na would be replaced by the width of 
the plate (perpendicular to the direction of propagation) and appropriately 
subdivided between the element termination nodes (two termination nodes for 
actual plate elements and one termination node for beam elements with a plane 
strain correction on the Young's modulus). 


Submerged Cylindrical Shell 

• determination of edge loading moments and forces 

The appropriate moments and forces employed to drive the shell edge employ 
the same equation (23) formulas, except that the propagating wave number root , y 
must be determined from equations (12b) in this case. 

• determination of shell edge termination dampers 

Again equations (27) can be used, except that y ss determined from equa- 
tions (27) is employed. 
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The appropriate fluid loading can be obtained by substituting equation 
(5a) into (5) (and converting into the -iyz + wt form by letting y = -y and 
<*) = -w) to obtain 


p(r,z) 


-p Wn to 2 H&j (rWk 2 -y 2 ) g-iyz + iwt 

k 2 -y 2 Hp (a^T k 2 -y 2 ) 


(29) 


Employing the large argument Bessel function approximations (ref. 5) of 


H^(X) s 'JttThxT e 1 (X _TT V 4 ) X ) s Hq ( X ) 

the equation (29) can be . ewritten for, \r*J y2 - kF | 

3 -(r/a-l)(a Nfy 2 -k 2 


p(r.z) z^lS^LI 


-iyz 


•sTp-k 1 ”' •vfr7<r 


e-^ 2 

> 10.0, in the form 

+iojt 

e 


(30) 


(31) 


where equation (31) can be used as the enforced pressure at the face z = 0. 

Equation (31) can also be used to drive plate models by first making a 
change in variable r = x + a (where x is the outward distance measured from the 
middle surface of the shell) and then appropriately passing to the limit as a+°° 
to obtain 


P(x,z) 


-Wq pu) 2 e 


-x 




e -iyz e +iwt 


(32) 


• determination of fluid loading on face z = L 

On this face we choose to absorb the incident pressure wave with an 
absorber type boundary condition rather than load it explicitly with equation 
(29) evaluated at z = L. This approach is consistent with the manner in which 
the shell was terminated. The relationship between pressure and applied forces 
in a finite element formulation is that the normal gradient of the pressure is 
proportional to the finite element model force. Therefore, the gradient of 
pressure normal to the cut, z = L is needed and given by (with equation (5)i(+yz-wt) 
arguments) 

= iyp(r.z) (34) 

and sinc° for steady state problems. 


ft 

it follows 


= n = rX A 

oZ co Bt a) y 


( 35 ) 
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is the boundary condition to employ along the vertical fluid cut z = L. For 
direct acoustic fluid elements (ref. 6), the finite element nodal force is 
related to the normal gradient by the relation 

F n ' i"™ 1 & ds (36) 

where [N] is the element shape function. For lumped force applications, equa- 
tion (36) is represented as, (where n‘ is the unit normal to the fluid and A A 
is an appropriate amount of area surrounding the node at the point of the F n 
appl i cation ) 


F 

n 



(37) 


and upon substituting equation (35) into (37) one obtains 


-AA- P 

U) 


(38) 


and represents the nodal force on the fluid as being proportional to the pres- 
sure time rate of change. Thus the internal force in a damper attached to the 
node would be equal and opposite in value and given by 


damp 

^n 



P 


where the damper value would be 


f-P = AAy direct pressure method damper 

u d co (no analogy) for z = L face 


(39) 


In the case of NASTRAN, the acoustic fluid is modeled by way of an analogy 
which reduces the elasticity equations to the Helmholtz equation by way of 
introducing dummy constants into the material matrix and allowing the displace- 
ment variable, u z , occupy the pressure variable (see ref. 3 for details). The 
analogy has been adopted herein, consequently, an additional factor of Pc2 
appears in equation (37), and consequently the frequency dependent NASTRAN 
axial damper needs to be modified by this same factor and thusly equation (39) 

is replaced by (replacing u >-to. y -*-y converts to NASTRAN convention; this results in 
the same damper si nee signs cancel ) 


C j = 


_ A MV 

w 


NASTRAN (ref. 3) analogy damper , . 

p for z = L face, with y directly from 

appropriate characteristic equation. 

It is noted, that so long as yis a real quantity (e.g., characteristic root 
equation for Y is obtained from equation (12b)), the termination condition in 
the z direction is a pure damper. However, where Yis complex (e.g., root from 
equation (12a)) substitution of y into the equation (40) damper will result in 
a complex damper. The real part can be treated as above, but the complex part 
will end up looking like a resistance (spring) when combined with the +ico 
appearing as a multiplier in the assembled damping matrix in NASTRAN. Thus, 
part of the damper (if it is present) can be applied as a nodal 
spring constant is 


the imaginar 
spring whose 


\ (-w) 


■) 


pc 2 j(+iw) = -aA pc 2 y 1 


NASTRAN. analogy spring 
(with y^ directly from 
appropriate characteristic 
192 equation) 
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wherein for decaying waves, y 1 will be a positive number making the spring 
constant positive. 

Equation (40) can equally be applied to plate mode) problems, upon properly 
interpreting aA. 

• determination of fluid loading on outer face r = r Q 

Here on the outer radial fluid face, the boundary condition is developed just 
as in the previous case, except that now the gradient has to be taken in the 
direction of the outward normal (n' = r). Thus operating on equation (29) with 
a( ) /a r and then employing the equation (30) large argument relations, the 
following result is obtained 

|jr = -*J7 2 -k 2 ^p(r,z) (41) 

Unlike equation (34), this result does not have a factor i rel a ting the gradi- 
ent of p and p itself, consequently, upon substition of eqirauon (41) into 
equation (37), one ob tains 

F n = -AA s[y T ^k I ~' p (42) 

and represents the nodal force on the fluid as being proportional to the pres- 
sure. Thus the internal force in a spring attached to the node would be equal 
and opposite in value and is given by 

= +AA vTy^-k 2 1 p 

where the spring constant value would be 


K p = aA •nTy 2 -k 2 


direct pressure method spring 
(no analogy) for r = ij, face 


(43) 


and in the case of employing the ref. 3 pressure analogy, the counterpart of 
equation (43) is 

K P = aA pcW7 z -k 2 ' NASTRAN (ref. 3) analogy (44) 

spring for r = r Q face 

When y is a root of equation (12b) (y% > k2), equation (44) implies a real 
spring is all that is needed. However, when y2 < k2 , and equation (12a) is 
employed, the roots, y, are, in general, complex which implies that equation 
(44) will have a real and imaginary part. The real part of KP is still a 
spring, however the imaginary portion can be converted into a damper by divid- 
ing the imaginary part coefficient (not including i) of equation (44) by u) (to 
compensate for the u multiplied in by the complex stiffness formation). 


Equation (44) can equally be applied to plate model problems, upon appropri 
ately interpreting aA. 
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Limitations of Membrane and Flexural Wave Modeling 


The scope of this study is limited to a driving frequency region where the 
influence of plate rotary inertia and shear deformation on the flexural motions 
are not important. In fact, this is the limitation of the propagating wave- 
length formulas given earlier. The ability of NASTRAN to account for these 
higher order effects is a topic all in itself and will not be addressed in this 
paper wherein stiffness and mass formulations along the lines given in ref. 7 
must be considered. 


Although we are interested in both plate and shell propagating waves, the 
classical plate theory (CPT) can provide insight with regard to the point at 
which theabove mentioned correction factors are needed. The CPT is defined 
here as the special case of equation (1) with a-* ®. Midlin’s classical paper 
(ref. 4) provides a guideline for when the characteristic wave equations pre- 
sented herein no longer predicts the proper waves number roots, y. As pointed 
out in ref. 4, the exact elasticity solution to the same problem shows that CPT 
can predict the correct flexural wave numbers when the corresponding 2 tt / y 
wavelengths are long in comparison with the thickness of the plate. As the 
wavelength gets smaller, the exact elasticity solution traveling wave phase 
velocity (w/y), has as its upper limit, the velocity of the Rayleigh surface 
wave. Hence the CPT cannot be expected to give good results for the wave num- 
bers, y, as the driving frequencies get increasingly large. Reference 4 pro- 
vides a plot of propagating phase velocity Cy = w/y (non-dimensional ized with 
respect to the shear wave velocity c s = Gibs' ) versus the plate thickness 
(non-dimensional ized with respect to the flexural wavelength, 2-n/yfor a fixed 
Poisson's ratio (v= 0.5) and is reproduced here as Figure 3. The CPT traight 
line plot, curve II, is nothing more than equation (9) rewritten in the form 
c Y h 2ir 


c s 


y 


*^6(l-v) 


(44a) 


and is compared against curve I, the exact elasticity solution. Curves III, 
IV, and V are different combinations of adding rotary interior and shear cor- 
rection factors. Figure 3 suggests that the shear correction addition has the 
biggest impact on correcting the CPT case. The shear thickness parameter T2, 
in the NASTRAN CCONEAX elements appears to attempt to account for shear 
effects, however, this point has not been pursued with regard to attempting to 
make NASTRAN propagate flexural waves for h/Xy> 0.15 


The results in Figure 3 suggest that a h/Xy (plate thickness- to-flexural 
wavelength ratio) of limit of 0.15 be maintained, i.e., 

^<0.15 (45) 

Ay 

or equivalently, 

$S0,15 (45a) 

where y is the flexural wave num ber root of the problem at hand. The limit 
given by equation (45a) employed Figure 3 which is a plot based on a fixed 
Poisson's ratio (i.e., v = 0.5). The guideline formula can still be used for 
other values of Poisson's ratio due to the weak dependence of the plots on v . 

In fact, for the common case of v = 0.3, a value of h/Xy = 0.15 implies a 
Cy / c s ordinate of 0.4598 for the CPT case and a value or 0.4032 in the exact 
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theory case. However, for v * 0.5, a h/x y = 0.15 implies a c y /c s ordinate of 
0.5441 in the CPT case and a value of 0.4629 in the exact theory case (computed 
from equation (40) of ref. 4). Thus, illustrating the Figure 3 v = 0.5 case is 
an extreme case (e.g., 17.5% deviaton for = 0.5 compared to 14.0% for v = 
0.3) with regard to the point at which CPT deviates from the exact theory. 

Thus, to test inequality (45a), one simply employs the appropriate charac- 
teristic equation developed (e.g., equation (12b)) and after obtaining the 
flexural wave -oot Y, verify whether in equation (45a) is satisfied. 


Comments on Mass Matrix 

The issue of an appropriate mass matrix for dynamics problems has long 
been a topic of discussion in the literature on finite elements. Here we take 
the simpliest approach, namely that of a diagonal lumped mass matrix for the 
shell. Upon employing the lumped approach, for say CC0NAX elements, NASTRAN 
generates zero valued rotary inertia mass matrix contributions. To fill this 
void, the work by ref. 8 was employed wherein the suggested diagonal rotary 
inertia terms at a node would be given by 

V = TT" [j0 + (?)] {46 

for the single element contribution at a rotation degree- of- freedom node where 
m« is the total mass of the element, l is the element length, h the element 
thickness. Thus end point nodes with one element framing into the node employs 
equation (46) once and internal nodes with two elements framing into the same 
node apply s equation (46) twice. 


DEMONSTRATION PROBLEMS 

A series of limited, yet fully representative, demonstration problems are 
presented here to illustrate the use of the procedures discussed in the theo- 
retical section. First an in vacuo cylindrical shell is treated employing three 
different mesh distributions. Secondly, the same shell (employing the finer 
mesh distribution) is solved with external water present. All solutions are 
obtained with COSMIC SOL-8 employing the VAX computer version of NASTRAN. 


In Vacuo Cylinder Example 

An infinite cylindrical shell is excited and is propagating a constant 
amplitude harmonic wave in the plus z direction. The shell (shown in Figure 1) 
is cut at points A and B and is modeled as a 3.2 inch long, 20 element CC0NEAX 
finite element model shown in Figure 4 (in the in vacuo case considered here, 
the fluid is omitted in the model). The point A (incident side of the shell) 
is driven with equations (23) and is terminated with three boundary dampers 
sized according to equations (27). The problem parameters correspond to the 
following data 
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(47) 






a = 40.0” 
v = 0\3 


h = 0.20" , 

p s = 0.000048, 


E = 154,149.39 psi 
oj = 9749.99 rad/ sec 


wherein the corresponding wave number roots to the characteristic equation (20) 
are given by 


and 


y = 0.163949685 
y = 1.67707302 


for the membrane root 
for the flexural wave root 


(48) 


hence the membrane wavelength is 12.566 inches and the flexural wavelength is 
3.7465 inches. 

• flexural wave example 

The shell termination dampers at cut B are computed upon substitution of 
the equation (48) flexural wave root into equations (27) to obtain 


Cg = 4.8819909 


Cg = 13.7309605 


C N = 14.0269417 


(49) 


and the shell edge forces, for an arbitrarily selected driven amplitude of 
Wq = 1.0 x 10" 6 , at the driver cut end A are obtained by substituting the flex- 
ural wave root into equations (23) to obtain (supressing the exponential time 
factor where dampers are installed with DMIG cards ): 

-A 


M‘ z = 7.9827681 x 10-2 
N z = 6.175285 x 10-4 


QA = 0.1338768 e™/2 


(50) 


The rotary inertia mass entries are computed with equation (46). The corre- 
sponding NASTRAN input is given in Figure 5 for reference purposes. The exact 
solution is a constant amplitude displacement wave varying as 

Consequently, the exact solution is easy to plot and is simply a constant mag- 
nitude radial displacement of magnitude W 0 = 1 x 10-5 and whose phase angle 
rolls off linearly with z and reploting with spatial period 2tt/y. The exact 
solution versus the NASTRAN solution is given in Figure 6 and shows good 
agreement between the exact solution and the corresponding NASTRAN solution. 
Since the model element spacing is 0.16 Inches, the solution model Is a 23^ 
element per wavelength case. In order to investigate the influence of mesh 
size on solution accuracy, the same 3.2 inch model was run for two coarser 
meshes; one having 5 elements and the other having only 3 elements. The same 
dampers and edge loads are reapplied to the coarser models. The resulting 
solutions are shown in Figure 7 and as can be seen, the quality of the solution 
(i.e., the ability of the cylindrical shell to propagate flexural waves) has 
degraded over the finer mesh example, particularly for the three element shell. 
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FIGURE 5. TWENTY ELEMENT CYLINDRICAL SHELL EXAMPLE DATA FOR FLEXURAL WAVE RUN 



NON-DIMENSIONAL RADIAL SHELL DISPLACEMENT |w/W Q |, NON-DIMENSIONAL RADIAL SHELL DISPLACEMENT 
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(z/X Y ); AXIAL COORDINATE NON-DIMENSIONAL IZE BY FLEXURAL WAVE LENGTH 


FIGURE 6. RADIAL SHELL DISPLACEMENT VS. AXIAL COORDINATE (FINE MESH) 



(z/X y ); AXIAL COORDINATE NON-DIMENSIONAL IZED BY FLEXURAL WAVE LENGTH 


FIGURE 7. RADIAL SHELL DISPLACEMENT VS. AXIAL COORDINATE (COARSE MESH) 
FOR 6 AND FOR 3H FLEXURAL WAVE LENGTHS PER ELEMENT 
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• membrane wave example 
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Again, the 20 element model is considered, except in this case, we propa- 
gate a membrane wave along the +z direction of the same cylindrical shell con- 
sidered for the flexural wave example. The procedure is exactly the same here, 
except that the membrane root (first of equations (48)) is employed to compute 
the shell termination dampers and driving moments and forces. Thus, substitut- 
ing the membrane root into equations (27) and (23), the model dampers and edge 
loads are in this case: 

= 0.4772607 C n » 0.012828529 C M = 143.484479 
0Q N 

and 

= 7.6290549 x 10-4 qA = 1.2507811 X 10-4 e irr/2 
*N ^ = -29.909871 

The solution response for the dominant axial motion is plotted in Figure 8 and 
is non-dimensional ized with respect to the, U 0 , axial wave amplitude magnitude 
(by equation (15), |U 0 | = 21.37985 x 10-6 when W 0 = 1.0 x 10-6). The radial 
deflection, non-dimensional ized with respect to W 0 , is given in Figure 9. The 
quality of the solution is seen to be very accurate, as would be expected since 
the wavelength is substantially longer than the flexural wave example (e.g., 
(2tt/- 0163949)70. 16 = 239^ elements per wavelength. 


Submerged Cylinder Example (Flexural Wave) 

Again, the same 20 element cylindrical shell is considered, except here, 
the shell is submerged in water (external water only where c = 60,000. in/ sec 
and p= 0.000096). Substituting the equation (47) parameters into equation 
(7b), and solving for the wave number roots, it is found that 


Y = 0.164124528 
and 

Y = 2.51631787 

for the membrane root 
for the flexural root. 

(51) 

Comparing the second of equations (51) to the in vacuo wave number, it is noted 
that the presence of the fluid shortened the wavelength by a factor of z 2/3. 
Substituting the flexural root into equations (27) and (23), the model dampers 
and edge loads (for a W 0 = 1 x 10-6 shell amplitude) are given by 

C e « 7.32505101 

, Cq = 46.3811654 , C N = 9.34867391 


and 



= 0.179713448 

, Q a = 0.452216160 eiir/2 , 



= 2.7283602 x 10-4 


In addition, the z = L boundary cut damper is given through equation (40), the 
r = rg = 43.2 inch fluid boundary spring by equation (44) and the z - 0 pres- 
sure loading by equation (31). The pressure values are enforced by applying 
the stiff spring approach, where the stiff spring constant is sized to be 1000 
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(z/X Y ), AXIAL COORDINATE NON-DIMENSIONALIZED BY MEMBRANE WAVE LENGTH 
FIGURE 8. AXIAL SHELL DISPLACEMENT VS. AXIAL COORDINATE (FINE MESH) 



U/Xy), AXIAL COORDINATE NON-DIMENSIONALIZED BY MEMBRANE WAVE LENGTH 
FIGURE 9. RADIAL SHELL DISPLACEMENT VS. AXIAL COORDINATE (FINE MESH) 
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times the regular pressure element stiffnesses and the applied force is the 
pressure times the stiff spring. The actual data input is shown in Figure 10, 
where the obvious repetitive pattern to the numbering system of Figure 4, per- 
mits us to leave out most of the grid coordinate cards, element cards, and 
fluid-to-structure DMIG connection cards while leaving behind representative 
examples of each kind. The resulting response for the radial deflection is 
plotted in Figure 11 and the corresponding pressure in the fluid (at the sur- 
face of the structure) is plotted in Figure 12. Both the deformation and pres- 
sure are seen to track the exact solution closely. A graphical representation 
of the entire pressure field amplitude is shown In Figure 13 through employing 
the PATRAN fringe color plotting feature. Plotting data in this fashion shows 
bands of data having the same magnitude spread, as a single color. Narrow 
bands at the surface spreading outward radially to increasingly wider bands 
show the exponential type decay in the pressure field. 


CONCLUDING REMARKS 

The results in this paper demonstrate a procedure by which the NASTRAN 
computer program can be employed to check the ability of the NASTRAN program to 
model membrane and flexural waves existing in both in vacuo and submerged 
cylindrical shells and flat plates. The study is limited to a range of fre- 
quencies where rotary inertia and shear correction factors are not necessary to 
model the corresponding wave propagation. For the demonstration problems con- 
sidered, the wave propagation ability of the elements considered appears to 
fall off rapidly, once lessthan 6 elements per wavelength are considered. For 
example, in the 6 element per flexural wavelength problem, the worse nodal 
point magnitude was in error by 12 % for the radial deflection, whereas the 
error was 36.7% for the element per wavelength example. It is recommended 
that the user make his own test with regard to mesh fineness necessary to 
achieve a particular level of accuracy. For example, when the propagating wave 
root is in the neighborhood of a cutoff frequency (i.e., a condition where no 
propagating wave exists), finer meshes than experineced In the demonstration 
considered in this paper may be needed. For most cases experienced by the 
authors, however, 10 elements per wavelength appears to provide good results 
for properly modeling the wave propagation for flexure and membrane waves In 
the kinds of elements considered herein. 
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FIGURE 10. TWENTY ELEMENT SL3MERGED SHELL EXAMPLE DATA FOR FLEXURAL WAVE RUN 



|p/P 0 !, NON-DIMENSIONAL FLUID PRESSURE (at r = a) |w/W Q U NON-DIMENSIONAL RADIAL SHELL DISPLACEMENT 



ORIGINAL PAGE IS 



(z/Xy) ; AXIAL COORDINATE NON-DIMENSIONAL I ZED BY FLEXURAL WAVE LENGTH 

FIGURE 11. RADIAL DISPLACEMENT VS. AXIAL COORDINATE FOR 154 FLEXURAL 
WAVE LENGTHS PER ELEMENT (SUBMERGED SHELL) 



(z/Xy); AXIAL COORDINATE NON-DIMENSIONAL IZED BY FLEXURAL WAVE LENGTH 


FIGURE 12. SURFACE PRESSURE VS. AXIAL COORDINATE FOR 154 FLEXURAL 
WAVE LENGTHS PER ELEMENT (SUBMERGED SHELL) 
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